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The  existence  and  spatial  development  of  hydrocarbon  cool  flames  in  a  spherical  vessel  under  the 
influence  of  mass  and  thermal  diffusion  have  been  investigated  by  numerical  methods.  The  purpose  is  to 
examine  the  nature  of  the  interaction  of  the  physics  and  chemistry  that  may  drive  an  oscillatory  reaction. 
The  conditions  correspond  to  those  that  would  be  experienced  at  zero  gravity,  as  has  been  recently  put  to 
experimental  test.  Comparisons  and  contrasts  with  responses  under  perfectly  mixed  conditions  are  made. 
The  numerical  simulation  was  based  on  a  skeleton  thermokinetic  scheme,  derived  from  that  of  Yang  and 
Gray,  in  a  three-variable  model  representing  two  intermediate  species  and  reactant  temperature.  Dirichlet 
and  Neumann  boundary  conditions  could  be  variously  selected.  The  equations  were  cast  in  one  dimension 
(spherical  symmetry)  and  integrated  using  the  numerical  algorithm  group  routine  D03PSF.  The  reactor 
surface  was  assumed  to  be  inert. 

Both  sustained  oscillatory  (i.e.,  multiple)  and  damped  cool  flames  were  predicted  to  exist  under  spatially 
uniform  conditions  resembling  those  reported  in  previous  experimental  studies.  The  phase  relationship 
between  the  chemical  species  and  temperature  in  sustained  oscillation  is  demonstrated.  There  was  strong 
evidence  for  the  negative  temperature-dependent  features.  No  sustained  oscillations  were  predicted  to 
occur  under  the  effect  of  diffusive  fluxes,  although  highly  damped  oscillations  were  still  able  to  exist.  This 
is  compatible  with  the  initial  experimental  observations  in  microgravity  conditions .  The  spatial  development 
reveals  the  growth  and  decay  of  the  reactive  intermediate  concentrations,  with  a  corresponding  expansion 
of  a  combustion  front  from  the  center  of  the  reaction  system  to  the  edge.  High  concentrations  of  inter¬ 
mediates  were  sustained  in  the  cooler  periphery  where  reaction  continued  to  be  supported.  Only  at  ab¬ 
normally  high  mass  diffusive  fluxes  could  sustained  oscillatory  reaction  be  recovered.  The  dependence  of 
oscillations  on  the  magnitude  of  mass  and  thermal  diffusion  coefficients  is  explored. 


Introduction 

Oscillatory  (or  multiple)  cool  flames  of  hydrocar¬ 
bons  and  other  organic  materials  in  oxygen  can  be 
generated  in  closed  or  flowing  systems  at  subatmos- 
pheric  pressures,  generally  at  vessel  temperatures  in 
the  range  500-700  K,  according  to  the  reactant.  The 
accompanying  temperature  changes  maybe  as  much 
as  200  K.  The  lowest  possible  reactant  pressures  for 
their  occurrence  are  not  less  than  about  20  kPa  (~  1/ 
5  atm).  This  means  that  under  terrestrial  conditions, 
self-heating  generates  an  accompanying  natural  con¬ 
vection  as  a  result  of  gravitational  effects.  Some  years 
ago,  Griffiths  et  al.  [1]  showed  that  the  buoyancy 
produced  complications  in  the  qualitative  structure 
of  the  observations  and  rendered  quantitative,  theo¬ 
retical  interpretation  of  the  nonlinear,  thermokinetic 


interactions  extremely  difficult,  if  not  impossible. 
They  also  advocated  the  use  of,  and  developed  with 
others,  well-stirred  systems  both  in  closed  and  flow¬ 
ing  form  [1,2],  Such  systems  have  proved  theirworth 
to  the  fundamental  understanding  of  the  oscillatory 
phenomena,  especially  because  numerical  model¬ 
ing — at  the  interface  between  theory  and  experi¬ 
ment — is  most  readily  applied  under  spatially  uni¬ 
form  conditions.  Indeed,  the  spatially  uniform 
criterion  is  a  prerequisite  for  the  simulation  of  the 
non-isothermal  oxidation  of  organic  compounds 
when  detailed  kinetic  models  involving  enormous 
numbers  of  species  are  exploited  [3]. 

The  heat  transport  process  under  spatially  uniform 
conditions  involves  transfer  at  the  interface  between 
the  gaseous  reactants  and  the  vessel  surface  and  is 
characterized  theoretically  and  numerically  by  a 
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(Newtonian)  heat  transfer  coefficient.  The  simplest 
foundations  for  thermal  ignition  theory  were  origi¬ 
nally  derived  by  Semenov  [4]  on  this  same  basis.  The 
alternative  criterion  for  heat  transport  in  thermal  ig¬ 
nition  systems,  that  of  pure  conductive  heat  loss,  was 
developed  quantitatively  by  Frank- Kamenetskii  [5], 
and  this  foundation  was  eventually  put  to  quantita¬ 
tive  experimental  test  by  Ashmore  et  al.  [6]  and  Fine 
et  al.  [7],  Such  tests  are  possible  in  gaseous  systems 
because  thermal  ignition  of  vapors  can  be  brought 
about  at  exceedingly  low  reactant  densities,  at  which 
the  Rayleigh  number  is  sufficiently  low  for  natural 
convection  to  be  suppressed. 

There  is  no  corresponding  opportunity  to  study 
oscillatory  cool  flame  phenomena  with  purely  con¬ 
ductive  heat  loss  under  terrestrial  conditions  and  so, 
perhaps,  there  has  been  little  impetus  to  investigate 
their  existence  and  associated  spatial  structure  either 
theoretically  or  by  numerical  analysis,  when  only  dif¬ 
fusive  fluxes  occur.  This  contrasts  with  the  lively  re¬ 
search  activity  that  prevails  with  regard  to  spatial 
structures  of  nonlinear,  isothermal,  kinetic  systems 

[8] ,  However,  an  impetus  has  emerged  recently  from 
experimental  studies  of  butane  +  oxygen  cool 
flames  under  microgravity  conditions  by  Pearlman 
and  co-workers  [9],  The  initial  studies,  on  butane  + 
oxygen  mixtures  in  a  closed  vessel,  include  an  image- 
intensified  photographic  record  of  the  cool  flame  de¬ 
velopment.  Although  not  unexpected,  a  striking  fea¬ 
ture  is  the  perfectly  symmetrical  development  of  the 
chemiluminescence  (CH20*)  from  the  center  of  the 
reaction  vessel,  as  an  expanding  circle  in  the  two- 
dimensional  field  of  view.  ( .1 120*  is  formed  only  in 
bimolecular  radical  interactions,  so  its  intensity  is  an 
indicator  of  where  high  radical  concentrations  exist. 
This  development  of  luminescence  contrasts  with 
the  buoyancy-driven  planar  front  reported  by  Grif¬ 
fiths  et  al.  [1]  and  detected  also  by  Pearlman  et  al. 

[9]  in  their  own  terrestrial  studies  which  comple¬ 
mented  the  near  zero-gravity  experiments.  Pearlman 
and  co-workers  also  found  in  the  initial  study  that 
whereas  multiple  (or  oscillatory)  cool  flames  propa¬ 
gated  under  terrestrial  conditions,  only  one  cool 
flame  occurred  when  heat  loss  was  dominated  by 
thermal  diffusion. 

These  novel  experiments  inspired  us  to  explore, 
by  numerical  analysis,  the  nature  of  the  interactions 
involved  in  thermokinetic  oscillations  when  the 
chemistry  is  coupled  both  to  heat  and  mass  transport 
by  diffusive  fluxes.  To  focus  on  the  physical  inter¬ 
actions,  in  particular  to  examine  the  way  in  which 
diffusive  fluxes  may  control  oscillatory  events,  we 
have  limited  the  chemistry  to  a  skeleton  kinetic 
scheme.  The  primary  interest  under  non-isothermal, 
but  non-adiabatic,  conditions  resides  in  how  a  neg¬ 
ative  temperature  coefficient  of  reaction  rate  con¬ 
trols  the  existence  of  oscillatory  reaction,  and  the 
way  in  which  the  oscillatory  behavior  may  be  af¬ 
fected  by  the  mode  and  magnitude  of  the  heat  and 


mass  transport  processes  involved.  To  this  end,  in 
the  present  paper  we  address: 

1.  The  behavior  of  a  simple,  kinetic  model  for  spa¬ 
tially  uniform,  non-isothermal  conditions. 

2.  How,  using  a  one-dimensional,  numerical  code, 
the  same  model  behaves  when  it  is  exposed  to 
mass  and  thermal  diffusion  in  spherical  symme- 
tiy. 

3.  The  spatial  structure  that  develops  under  the 
conditions  of  item  no.  2  at  Lewis  number  ( Le )  = 
1,  and  also  when  Le  is  varied. 

Kinetic  Foundation 

The  skeleton  thermokinetic  scheme  is  derived 
from  the  Yang  and  Gray  model  [10],  but  incorporat¬ 
ing  two  intermediate  species,  which  in  the  context 
of  alkane  oxidation  may  be  likened  to  a  free  radical 
propagator  (A)  and  a  peroxydic,  molecular  branching 
intermediate  ( B ).  There  are  competitive  branching 
and  termination  reactions  with  different  activation 
energies.  All  of  the  heat  is  regarded  to  be  generated 
at  the  branching  cycle.  No  supplementary  propaga¬ 
tion  processes  are  considered,  which  might  be  the 
main  source  of  heat  release  in  reality.  P  is  a  reactant 
(or  precursor),  and  C  represents  final  products.  One 
distinction  of  the  role  of  the  precursor  from  that  of 
a  hydrocarbon  fuel  is  that  P  is  deemed  not  to  be 
involved  in  any  subsequent  steps  [10,11],  This  means 
that  the  initiation  step  can  be  slowed  sufficiently  that 
the  computed,  time-evolved  states  approximate  well 
to  stationary  states  of  the  system  because  there  is  no 
significant  perturbation  by  marked  reactant  con¬ 
sumption.  Damped  oscillatory  modes  are  readily  dis¬ 
tinguished  from  sustained  oscillatory  states,  in  con¬ 
sequence.  No  activation  energy  need  be  assigned  to 
the  initiation  step  when  this  approach  is  invoked. 
The  negative  temperature  dependence  of  reaction 
rate  (ntc),  which  is  required  in  a  simulation  of  alkane 
oxidation,  is  conferred  by  the  competition  of  a  ter¬ 
mination  reaction  which  has  a  higher  activation  en¬ 
ergy  than  the  branching  reaction  (Etl  >  Eb).  An  aged 
silica  vessel,  as  used  in  the  microgravity  study,  is  ex¬ 
pected  to  be  relatively  inert  to  surface  termination 
of  intermediates.  A  weak  termination  of  species  B 
was  included  to  represent  this  state,  the  magnitude 
of  which  (At2)  determines  the  lower  temperature 
bound  for  oscillations.  Arrhenius  parameters  were 
chosen  for  the  data  set  to  tune  the  reactive  region 
to  be  comparable  with  that  of  experimental  obser¬ 
vations,  and  the  exothermicity  was  chosen  to  signify 
the  modest  heat  output  associated  with  low  tem¬ 
perature  oxidation  chemistiy.  The  kinetic  scheme 
and  its  parameters  are  shown  in  Table  1. 

This  scheme  has  attributes  of  the  isothermal  cubic 
autocatalytic  model  developed  by  Gray  and  Scott 
[11],  the  autocatalytic  step  of  which  is  replaced  by  a 
highly  nonlinear,  non-isothermal  branching  step.  An 
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TABLE  1 

Kinetic  scheme  and  parameters 


Pre-exponential 
Factor  (s-1) 

Activation  Energy 
(kj  mol-1) 

Exothermicity 
(Q/J  mol-1) 

Initiation 

P  — >  A 

Aj  =  1  X  10  6 

Ei  =  0 

0 

Propagation 

A  — >  B 

Ap  =  1  X  106 

Ep  =  0 

0 

Termination  1 

A  -»  C 

Atl  =  1  X  1011 

£.i  =  75 

0 

Branching 

B  — >  2A 

Ab  =  8  X  103 

cr 

II 

Or 

1  X  10s 

Termination  2 

B  — >  C 

At 2  =  0.60 

Ea  =  0 

0 

energy  conservation  equation  is  included  here  in  ad¬ 
dition,  yielding  a  three-variable  model,  represented 
by  the  concentrations  of  A  and  B  and  the  reactant 
temperature  T.  The  relevant,  fully  dimensional  con¬ 
servation  equations  are  as  follows,  where  the  ele¬ 
mentary  rate  constants  are  represented  by  k  and  the 
concentrations  of  the  respective  species  are  denoted 
by  lowercase  letters. 

^  =  kp„e{~klt)  -  ktl(T)a  -  k?a 

+  2  kh(T)b  +  DAV2a  (1) 

-  =  V  -  kb(T)b  -  kt3b  +  DBV2b  (2) 

at  r 

cdT  Qkh(T)b  1  (T  -  Ta)  „„ 

-  =  v  h3  -  -2 - +  DtV2T  (3) 

df  Cva  fN 

The  Lewis  number,  Le,  is  defined  as  DT/DA.  Both 

Newtonian  and  diffusive  thermal  heat  transport 
terms  are  included  in  the  energy  conservation  equa¬ 
tion.  These  were  flagged  to  permit  adoption  of  either 
spatially  uniform  or  diffusive  heat  transport  condi¬ 
tions  in  the  numerical  procedure  described  in  the 
next  section.  An  additional  option  in  the  diffusive 
flux-driven  calculations  was  to  choose  boundary  con¬ 
ditions  that  reflected  either  total  loss  or  no  loss  of 
the  autocatalytic  intermediate  at  the  surface.  The 
boundary  conditions  were,  in  each  case, 

da  3T 

—  =  0  and  —  =  0  at  x  =  0  4 

dx  dx 

dT 

—  =  0  at  x  =  r, 
dx 

for  spatially  uniform  conditions  (5) 
(T  —  Ta)  =  0  and  (optionally)fc  =  0  at  x  =  r, 
for  purely  conductive  conditions  (6) 

Heat  capacity  and  reactant  density  were  taken  to 
be  105  J  mol-1  K  1  and  8.2  X  10“ 6  mol  cm'3, 
respectively,  and  the  (constant)  volumetric  heat  ca¬ 
pacity  was  8.6  X  10“4  J  m“3  K_1.  The  chosen  den¬ 
sity  corresponds  to  an  initial  reactant  pressure  of  40 


kPa  at  600  K,  which  is  typical  of  the  experimental 
conditions  [9],  The  numerical  results  are  not 
strongly  sensitive  to  the  heat  capacity  and  not  at  all 
sensitive  to  the  reactant  density  because  it  is  not  in¬ 
volved  in  subsequent  steps.  The  requirement  for  ox¬ 
ygen  may  be  regarded  to  be  subsumed  in  the 
(pseudo)  first-order  rate  parameters.  The  simula¬ 
tions  represent  experiments  in  a  500  cm3  spherical 
vessel. 

Numerical  Methods 

Equations  1-3,  subject  to  boundary  conditions  4- 
6,  were  cast  in  non-dimensional  form,  representing 
the  classical  one-dimensional  shapes  (infinite  slab, 
infinite  cylinder,  or  sphere),  and  integrated  using  the 
numerical  algorithm  group  (NAG)  routine  D03PSF 

[12] ,  This  routine  is  set  up  for  the  integration  of  a 
system  of  nonlinear  convection-diffusion  equations 
in  one  dimension,  using  the  method  of  lines  to  re¬ 
duce  the  partial  differential  equations  to  ordinary 
differential  equations  (ODEs).  The  ODE  system  is 
then  solved  by  a  backward  difference  (BDF)  method 

[13] ,  The  BDF  algorithm  makes  use  of  formulae,  the 
order  of  which  are  varied  automatically  up  to  order 
5.  The  scheme  is  an  implicit,  variable  step-size 
method,  the  time  step  and  order  being  selected  au¬ 
tomatically  to  satisfy  a  user-specified  tolerance.  This 
allows  temperature  spikes,  which  may  develop  on 
very  short  time  scales,  to  be  accurately  modeled 
while  not  restricting  the  time  step  in  other  phases  of 
the  simulation.  A  modified  Newton  iteration  scheme 
is  used  to  solve  the  system  of  nonlinear  equations  at 
each  time  step. 

The  problem  was  discretized  using  a  finite  volume 
approach  [12].  The  source  terms  were  calculated  at 
each  node,  and  the  diffusive  terms  were  evaluated 
at  each  nodal  midpoint  using  values  derived  from  the 
nodes  at  either  side  [14],  These  values  were  calcu¬ 
lated  by  assuming  that  the  solution  varied  within  a 
cell  and  employed  standard  upwind  techniques  com¬ 
bined  with  a  slope  limiter  [15],  In  the  present  ap¬ 
plication,  the  domain  was  meshed  uniformly  with 
126  points.  Neumann  conditions  were  specified  at 
the  symmetry  boundary,  and  a  mixture  of  Neumann 
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Fig.  1.  Temperature-time  profiles  from  the  integration 
of  equations  1-3  under  spatially  uniform  conditions  at  fN 
=  0.2  s.  The  initial  and  boundary  temperatures  are  re¬ 
spectively  650-750  K  in  20  K  increments.  This  spans  the 
transition  from  sustained  oscillations  through  damped  os¬ 
cillatory  reaction  to  stable  stationary  state  reaction.  The 
sustained  oscillations  set  in  at  Ta  =  590  K. 


0.0-1 - , - 1 - . - 1 - , - 1 - . - 1 - . - 1 

650  700  750  800  850  900 
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Fig.  2.  A  phase  plane  map  representing  the  relationship 
between  the  concentration  of  intermediate  B  and  reactant 
temperature  over  the  time  interval  7. 8-9. 5  s  at  Ta  =  650 
K  (see  Fig.  1).  The  rise  and  fall  of  [B]  is  predicted  to  pre¬ 
cede  that  of  T,  as  marked  by  the  arrows. 

and  Dirichlet  conditions  were  specified  at  the  outer 
boundary,  as  noted  above. 

The  code,  D03PSF,  provides  accuracy  up  to  order 
5  for  a  single  time  step.  For  a  number  of  time  steps, 
the  local  error  tolerances  were  investigated,  and  val¬ 
ues  were  chosen  that  did  not  compromise  the  solu¬ 
tion.  For  spatial  accuracy,  the  solutions  from  differ¬ 
ent  numbers  of  mesh  points  were  compared.  The 
chosen  number  reflected  a  satisfactory  compromise 
between  computational  speed  and  solution  accuracy. 

Results  and  Discussion 

Spatially  Uniform  Conditions 

A  Newtonian  cooling  time  (fN)  of  0.2  s  was  as¬ 
signed  to  the  system,  which  is  a  reasonable  repre¬ 
sentation  for  a  well-stirred  closed  system  [16],  The 


results  of  the  simulations  for  temperature  versus 
time  over  a  vessel  temperature  range  650-750  K  are 
shown  in  Fig.  1.  The  onset  of  sustained  oscillations 
occurs  at  590  K  (not  shown),  which  is  set  by  the 
magnitude  of  At2,  but  there  is  a  transition  to  a 
damped  oscillatory  state  at  about  Ta  =  660  K  (Fig. 
1).  The  damped  oscillations  converge  to  a  fixed  sta¬ 
tionary  state  temperature  (T  =  782  K),  correspond¬ 
ing  to  that  attained  in  the  stable,  non-oscillatory 
high-temperature  state,  at  Ta  >  740  K.  The  fixed 
value  for  the  final  reactant  temperature  is  consistent 
with  the  response  to  an  ntc  of  reaction  rate  in  hy¬ 
drocarbon  combustion,  as  observed  experimentally 
in  the  combustion  of  alkanes  following  rapid  com¬ 
pression  [17]  and  in  propane  combustion  in  a  high- 
pressure  flow  reactor  [18]. 

Phase  relationships  can  be  distinguished  when  sta¬ 
ble  oscillations  occur,  as  shown  in  the  two-dimen¬ 
sional  phase  plane  for  the  relationship  between  the 
concentration  of  B  and  the  temperature  T  during 
oscillations  at  650  K  (Fig.  2).  The  maximum  in  the 
concentration  of  B  precedes  the  maximum  of  T  in 
each  cycle.  The  concentrations  of  the  intermediates 
A  and  B  rise  and  fall  in  phase,  and  so  are  not  shown. 

Spatially  Non-uniform  Conditions 

Spatially  non-uniform  behavior  was  generated  as 
a  result  of  diffusive  fluxes  for  mass  and  energy,  taking 
Da  =  Db  =  3  cm2  s-1  and  Le  =  1,  with  no  loss  of 
B  following  diffusion  to  the  surface.  The  diffusion 
coefficient  would  not  normally  be  expected  to  ex¬ 
ceed  this  value  [5],  Sustained  oscillations  were  not 
detected  at  any  temperature  under  these  conditions, 
the  only  manifestation  being  that  of  highly  damped 
oscillations  to  a  stationary  state  (Fig.  3).  The  tem¬ 
perature  records  in  Fig.  3  relate  to  behavior  at  the 
center  of  the  vessel.  Further  study  of  the  domain  of 
Da  showed  that  sustained  oscillations  could  be  es¬ 
tablished  under  the  effect  of  diffusive  fluxes,  for  ex¬ 
ample,  at  the  hypothetical  case  at  DA  >  15  cm2  s  “ 1 
(Le  =  1)  when  Ta  =  590  K,  but  rising  to  DA  >  30 
cm2  s_1  at  Ta  =  650  K  (see  later).  An  experimental 
pressure  record  from  the  microgravity  experiments 
is  also  shown  in  Fig.  3.  There  is  a  relatively  slow 
entry  of  reactants,  during  which  time  the  reaction 
develops  to  a  damped  oscillatory  state,  exemplified 
by  the  initial  overshoot  in  pressure,  followed  by  a 
limited  number  of  smaller  amplitude  oscillations  [9] . 

The  spatial  variations  associated  with  the  devel¬ 
opment  of  reaction  were  also  investigated  with  re¬ 
spect  to  temperature  and  intermediate  product  con¬ 
centrations  over  a  range  of  conditions.  These  are 
illustrated  in  Figs.  4  and  5  at  the  threshold  of  the 
sustained  oscillatory  mode,  for  which  Ta  =  590  K  at 
Le  =  1 .  The  response  of  the  three  variables  over  the 
time  interval  12.52-13.41  s,  corresponding  to  the 
initial  temperature  excursion,  is  shown  in  Fig.  4.  The 
temperature  develops  earliest  at  the  center  of  the 
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Fig.  3.  Temperature-time  profiles  from  the  integration 
of  equations  1-3  in  one  dimension  (spherical  symmetry) 
under  diffusive  heat  and  mass  fluxes,  DA  =  DB  =  DT  = 
3  cm2  s_1.  The  initial  and  boundary  temperatures  are,  re¬ 
spectively,  590-750  K  in  40  K  increments.  This  spans  the 
full  range  from  the  onset  of  damped  oscillatory  reaction  to 
stable  stationary  state  reaction.  No  sustained  oscillations 
were  predicted  at  these  conditions.  An  experimental  pres¬ 
sure  record  from  the  microgravity  experiments  is  also 
shown.  This  was  obtained  from  the  mixture  0.25  C4H10  + 
0.25  02  +  0.5  Ar  entering  a  reaction  vessel  at  583  K  to  an 
initial  pressure  of  29  kPa. 


sphere,  but  the  successive  temperature  profiles 
(12.85-13.27  s)  convey  a  sense  of  the  outward  de¬ 
velopment  of  a  combustion  front  (Fig.  4a).  This  is 
not  a  self-sustained  combustion  wave  insofar  that  it 
does  not  propagate  fully  to  the  surface. 

Both  A  and  B  exhibit  maxima  in  their  concentra¬ 
tions  which  are  displaced  outward  throughout  the 
time  period  of  the  thermal  records  in  Fig.  4a  (Fig. 
4b  and  c).  The  maxima  of  A  and  B  are  not  coincident 
in  time  and  space,  however.  The  progress  of  the  peak 
concentration  of  A  toward  the  edge  of  the  reaction 
vessel  (Fig.  4b)  is  reminiscent  of  the  symmetrical 
expansion  of  the  chemiluminescent  region  observed 
experimentally  [9] .  The  intermediates  are  almost  en¬ 
tirely  consumed  at  the  center  of  the  system  in  the 
later  stages.  Whereas  B  is  able  to  accumulate  at  the 
surface  to  yield  its  highest  concentrations  there  in 
the  late  stages  of  the  thermal  cycle,  a  marked  gra¬ 
dient  in  the  concentration  of  A  is  predicted  at  the 
surface.  The  minor  oscillations  in  Fig.  4b  are  a  nu¬ 
merical  artifact. 

The  evolution  during  the  time  period  from  13.69 
s  to  14.81  s  is  depicted  in  Fig.  5a-c,  and  this  corre¬ 
sponds  to  the  interval  during  which  the  center  tem¬ 
perature  relaxes  from  its  maximum  at  about  1000  K 
to  around  890  K  (see  Fig.  3).  Just  at  the  end  of  this 
time  period,  there  is  a  thermal  response  to  chemical 
activity  at  about  2.5  cm  along  the  radius  (Fig.  5a). 
This  distorts  the  thermal  profile  to  a  virtually  flat 
central  region  (f  =  15.09  s)  which,  in  the  subsequent 
stage  of  reaction,  evolves  to  a  slight  dip  in  tempera¬ 
ture  toward  the  vessel  center  ( t  =  15.51  s).  There 
is  no  substantial  qualitative  change  throughout  the 


Fig.  4.  The  predicted  evolution  of  the  spatial  distribu¬ 
tions  of  (a)  temperature  T,  (b)  intermediate  species  A,  and 
(c)  intermediate  species  B  across  a  radius  of  the  spherical 
mass  at  selected  times  during  reaction  at  Ta  =  590  K.  The 
respective  times  are  marked,  and  they  correspond  to  the 
period  for  the  development  of  the  initial  temperature  peak. 
The  minor  fluctuations  in  the  profiles  for  A  arise  from  the 
choice  of  output  frequency  of  data  points  and  the  mesh. 
The  broken  fines  signify  profiles  for  which  there  has  been 
a  change  of  direction  in  the  temporal  evolution. 

respective  spatial  concentration  profiles  of  A  or  B 
(Fig.  5b  and  c),  and  the  growth  in  chemical  activity 
is  confined  to  the  outer  regions  of  the  reactor.  Once 
the  initial  thermal  gradient  has  developed,  mass  and 
thermal  diffusion  seem  not  to  be  sufficient  to  affect 
the  behavior  in  other  regions. 
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Fig.  5.  The  predicted  evolution  of  the  spatial  distribu¬ 
tions  of  (a)  temperature  T,  (b)  intermediate  species  A,  and 
(c)  intermediate  species  B  across  a  radius  of  the  spherical 
mass  at  selected  times  during  reaction  at  Ta  =  590  K.  The 
respective  times  are  marked.  They  correspond  to  the  pe¬ 
riod  during  which  there  is  a  relaxation  of  the  first  tem¬ 
perature  peak  and  development  of  the  second  one.  A 
slightly  later  temperature  profile  is  also  shown  in  which  the 
temperature  at  2.0  cm  exceeds  that  at  the  center.  The  bro¬ 
ken  lines  signify  profiles  for  which  there  has  been  a  change 
of  direction  in  the  temporal  evolution. 


Fig.  6.  Temperature-time  profiles  from  the  integration 
of  equations  1-3  in  one  dimension  (spherical  symmetry) 
under  the  respective  diffusive  heat  and  mass  fluxes.  Full 
line,  Da  =  Db  =  16  cm2  s  1 ,  DT  =  3  cm2  s  1  (Le  = 
0.1875);  broken  line,  DA  =  DB  =  DT  =  16  cm2  s  1  (Le 
=  1);  dotted  line,  DA  =  DB  =  3  cm2  s  f  DT  =  16  cm2 
s_1  (Le  =  5.33).  The  initial  and  boundary  temperatures 
are  590  K.  These  results  support  the  interpretation  that  the 
mass  diffusion  is  the  dominant  term  in  determining 
whether  or  not  sustained  cool  flame  oscillations  can  occur 
under  the  influence  of  diffusive  fluxes. 

To  what  extent  then  does  the  limitation  on  mass 
or  thermal  diffusion  cause  the  failure  of  sustained 
oscillations?  The  enhancement  of  both  diffusive  pro¬ 
cesses  provokes  sustained  oscillations,  as  shown  in 
the  thermal  record  at  Dx  =  16  cm2  s  ,  Le  —  1 
(Fig.  6).  There  is  the  same  qualitative  effect  when 
Dt  =  3  cm2  s-1  and  DA  —  16  cm2  s_1,  at  Le  = 
0.1875  (Fig.  6).  How  the  rate  of  heat  loss  by  the 
enhancement  of  Dx  affects  the  oscillatory  frequency 
can  also  be  seen  in  Fig.  6.  Reversing  the  diffusivities 
so  that  Dx  =  16  cm2  s_1  (Le  =  5.33)  allows  only 
highly  damped  oscillations  to  occur  throughout  the 
temperature  range.  Spatial  temperature  and  concen¬ 
tration  changes  for  B  associated  with  the  time  period 
13.70-15.30  s  at  Le  —  0.1875  are  shown  in  Fig.  7. 
Clearly,  in  these  circumstances  the  diffusion  of  B  is 
sufficiently  high  to  reduce  the  intermediate  species 
concentration  at  the  periphery  of  the  reaction  vessel, 
which  has  a  consequence  on  the  reaction  rate 
throughout  the  system. 

Conclusions 

When  a  perfectly  mixed  state  exists,  both  sus¬ 
tained  and  damped  oscillatory  modes  are  predicted 
to  occur  using  a  thermokinetic  model  to  represent 
hydrocarbon  oxidation  in  the  cool  flame  region.  By 
contrast,  heat  and  mass  transport  under  normal,  dif¬ 
fusive  conditions  are  not  sufficiently  high  to  promote 
sustained  oscillations.  The  diffusive  fluxes  have  to  be 
raised  to  rather  higher  values  for  sustained  oscilla¬ 
tions  to  be  induced.  Nevertheless,  damped  oscilla¬ 
tory  reaction  states  can  exist  at  normal  diffusivities, 
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Fig.  7.  The  predicted  evolution  of  the  spatial  distribu¬ 
tions  of  (a)  temperature  T  and  (b)  intermediate  species  B 
across  a  radius  of  the  spherical  mass  at  selected  times  dur¬ 
ing  reaction  at  Ta  =  590  K  for  the  conditions  DA  =  DB  = 
16cm2s”1,  Dy  =  3cm2s~1(Le  =  0.0625).  The  respective 
times  are  marked.  They  correspond  to  the  period  in  which 
the  second  cool  flame  develops  and  subsides.  The  broken 
lines  signify  profiles  for  which  there  has  been  a  change  of 
direction  in  the  temporal  evolution. 

which  are  a  manifestation  of  cool  flame  activity  and 
would  appear  visually  as  a  single  cool  flame. 

The  failure  to  generate  sustained  oscillations  is  at¬ 
tributed  to  the  accumulation  of  reactive  intermedi¬ 
ates  in  the  outer  reaches  of  the  reacting  mass,  where 
the  temperature  is  lowest  (Figs.  4  and  5),  so  that  the 
chemical  activity  is  then  sustained  in  these  localized 
regions  and  cooling  in  the  center  cannot  then  occur. 
The  dispersion  of  the  intermediates  by  enhanced 
mass  diffusion  causes  a  major  qualitative  change, 
such  that  sustained  oscillations  are  then  possible. 

An  alternative  route  to  reduction  of  the  interme¬ 
diate  concentration  in  the  periphery  of  the  vessel 
would  be  to  promote  surface  termination.  However, 
if  the  boundary  conditions  are  transformed  to  permit 
full  destruction  of  B  at  the  surface,  as  in  equation  6, 
this  so  diminishes  the  reactivity  that  strong  exother¬ 
mic  reaction  becomes  possible  in  reasonable  time- 
scales  only  at  Ta  >  650  K,  where  highly  damped 
oscillations  are  inevitable.  This  is  consistent  with  the 


experimentally  observed  suppression  of  cool  flames 
when  efficient  chain  terminating  surfaces  are  used 
[19].  There  was  also  no  indication  of  a  propagating 
front  when  these  boundary  conditions  were  applied. 

It  is  not  intended  to  imply  that  the  limiting  dif¬ 
fusion  coefficients  or  Lewis  numbers  given  here 
could  be  accessible  experimentally  or,  indeed,  that 
the  same  limits  would  necessarily  prevail  if  a  com¬ 
prehensive  kinetic  scheme  to  represent  normal  bu¬ 
tane  oxidation  were  investigated  numerically  in  a 
one-dimensional  model.  However,  we  believe  that 
these  are  the  first  indications  of  how  the  generation 
of  multiple  cool  flames  of  hydrocarbons  are  affected 
by  heat  and  mass  transport  as  a  result  of  diffusion. 
Confirmation  has  been  gained  in  subsequent  exper¬ 
iments  by  Pearlman  et  al.  when  helium  has  been 
used  as  a  diluent. 
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